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Abstract.  Deformation  induced  heating  of  reactive  solids  is  a  physically  complex 
process.  As  such,  the  effects  of  meso-structure,  component  thermomechanical  properties, 
component  mass  fractions,  and  porosity  on  their  impact  response  is  not  well-understood. 
In  this  study,  an  explicit,  2-D,  Lagrangian  finite  and  discrete  element  technique  is  used 
to  examine  thermomechanical  fields  in  metal-explosive  (aluminum-HMX)  particle  mix¬ 
tures  due  to  piston  supported  uniaxial  deformation  waves.  The  meso-scale  description  uses 
a  plane  strain,  thermoelastic-viscoplastic  and  friction  constitutive  theory  to  describe  the 
motion  and  deformation  of  individual  particles,  and  an  energy  consistent,  penalty  based 
method  to  describe  inter-particle  contact.  The  deformation  response  of  material  having 
an  initial  solid  volume  fraction  of  <j>°s  =  0.835  is  characterized  for  different  metal  mass 
fractions  and  wave  strengths.  Predictions  indicate  that  the  response  can  be  classified  into 
strength  dominated  and  pressure  dominated  regions  depending  on  wave  strength.  Average 
thermomechanical  fields  that  define  the  macro-scale  wave  structure  are  found  to  differ 
both  qualitatively  and  quantitatively  between  the  two  regions. 


1  Introduction 

Reactive  solids  often  consist  of  mixtures  of  high-explosive,  metal,  and  oxidizer  particles  (size  ~ 
1 0-200 yum)  embedded  within  a  polymeric  binder.  It  has  been  shown  that  inclusion  of  metal  in  ex¬ 
plosive  mixtures  can  affect  their  impact  sensitivity,  detonation,  and  post-detonation  behavior,  though 
it  remains  fundamentally  unclear  how  meso-structure,  component  (or  phase)  thermomechanical  prop¬ 
erties,  and  metal  mass  fraction,  affect  the  impact  sensitivity  of  the  mixture.  The  objectives  of  this 
study  are  two-fold:  1)  to  characterize  the  relative  importance  of  phase-specific  volumetric  and  surface 
dissipation  and  its  dependence  on  meso-structure,  metal  mass  fraction,  and  wave  strength;  and,  2) 
to  investigate  macro-scale  manifestations  of  the  particle-scale  response  by  spatially  averaging  meso- 
scale  fields  over  suitably  chosen  representative  volumes  (areas).  The  first  objective  is  important  for 
identifying  physical  mechanisms  responsible  for  ignition,  whereas  the  second  objective  is  important 
for  developing  macro-scale,  multiphase  constitutive  theories  to  accurately  assess  the  performance  of 
reactive  solids  in  practice  [1,2],  Multiphase  theories  are  routinely  posed  in  terms  of  phase  interaction 
terms  that  implicitly  account  for  particle-scale  phenomena  that  are  difficult  to  experimentally  resolve. 

The  initial  material  configuration  used  in  this  study  consists  of  circular,  randomly  packed  alu¬ 
minum  and  HMX  particles  having  an  average  diameter  of  d  «  60 gm,  as  illustrated  in  Fig.  1.  Inter¬ 
particle  voids  are  massless.  This  idealized  geometry  is  chosen  to  mitigate  singularities  that  would  occur 
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Fig.  1.  Illustration  of  an  aluminized  HMX  mixture  corresponding  to  A„,  =  0.2  and  Up  =  500  m/s.  The  blue 
particles  are  metal  and  the  red  particles  are  HMX. 

in  the  vicinity  of  sharp  particle  corners  that  are  impossible  to  numerically  resolve.  Nonetheless,  this 
simplified  initial  configuration  still  enables  leading-order  effects  of  viscoplastic  deformation  and  inter¬ 
particle  friction  to  be  assessed.  The  average  initial  solid  volume  fraction  of  the  material  is  <p°s  =  0.835. 
For  brevity,  only  predictions  for  piston  speeds  of  Up  =  100  and  500  m/s  are  discussed  as  they  enable 
differences  between  the  strength  and  pressure  dominated  deformation  regions  to  be  highlighted. 


2  Model  description 

The  model  problem  is  posed  as  a  multi-particle  contact  problem  for  the  evolution  of  thermomechanical 
fields  within  particles  and  interactions  between  particles.  Evolution  equations  for  the  momentum  and 
energy  of  each  particle  are  given  in  Lagrangian  form  by  the  following: 

pii  =  V  •  tr,  (1) 

pcj  =  -V  •  q  +  pr.  (2) 

Here,  p  is  the  local  density,  cr  is  the  Cauchy  stress  tensor,  r  is  the  deformation  induced  heating,  q  is  the 
heat  flux,  c„  is  the  constant  volume  specific  heat,  and  V  =  /-  *  is  the  spatial  gradient  operator.  Particles 
are  assumed  to  be  initially  stationary  and  stress  free.  To  satisfy  the  entropy  inequality,  a  Clausius 
Duhem  entropy  constraint  is  used  which  is  given  by: 

pr]  -  T~l(pr  -  V  q)  +  q  VT^1  >  0,  (3) 

where  r\  is  the  mass-specific  entropy.  Additionally,  the  following  contact  conditions  are  imposed  on 
particle  boundaries  T : 

cr  ■  n  =  tf  on  T  V  t,  (4) 

-kVT  ■  n  =  qf  +  qc  on  T  V  t.  (5) 

Here,  n  is  the  local  unit  normal  vector  to  T,  t,  is  the  local  contact  traction  vector,  qf  is  local  frictional 
heat  flux,  qc  is  the  local  heat  flux  needed  to  impose  ideal  thermal  contact,  and  k  is  the  thermal  con¬ 
ductivity.  The  stress-strain  behavior  of  the  particles  is  modeled  by  a  hyperelastic,  multiplicative,  finite 
strain  constitutive  model  with  a  Von-Mises  yield  criterion  and  viscoplastic  flow  rule,  and  friction  is 
modeled  using  an  Amo ntons -Coulomb  stick/slip  theory.  A  detailed  description  of  the  model  equations, 
contact  conditions,  and  constitutive  theory  is  given  in  Ref.  [3].  Key  thermomechanical  properties  for 
HMX  and  aluminum  used  in  the  simulations  are  given  in  Table  1 . 


2.1  Numerical  technique 

The  model  equations  are  computationally  solved  subject  to  the  boundary  conditions  using  a  com¬ 
bined  finite  and  discrete  element  technique.  The  finite-element  method  (FEM)  is  used  to  numerically 
integrate  the  unsteady  2-D  conservation  equations  and  the  viscoplastic  flow  rule  governing  inelastic 
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Table  1.  Thermomechanical  properties. 


Property 

HMX 

Aluminum 

Density  (kg/m3) 

1903 

2400 

Yield  Stress  (GPa) 

0.37 

0.25 

Poisson’s  Ratio  (-) 

0.20 

0.35 

Young’s  Modulus  (GPa) 

24.0 

70.0 

Bulk  Modulus  (GPa) 

13.33 

77.8 

Specific  Heat  (kJ/kg/K) 

1.50 

0.90 

Thermal  Conductivity  (W/m/K) 

0.5 

237.0 

Longitudinal  Sound  speed  (km/s) 

2.65 

5.69 

Shear  Sound  Speed  (km/s) 

2.29 

3.29 

deformation,  whereas  the  discrete  element  method  (DEM)  is  used  to  account  for  interactions  between 
particles.  The  DEM  utilizes  a  conservative  potential  based  penalty  method  in  which  the  normal  con¬ 
tact  traction  between  particles  is  estimated  by  penalizing  their  penetration,  and  frictional  tractions  are 
estimated  using  a  penalty  regularized  Amontons-Coulomb  law  [3].  The  initial  particle  configuration 
is  generated  using  a  pseudo-gravity  drop  method  that  consists  of  two  steps.  First,  an  initial  packing 
arrangement  is  generated  by  randomly  placing  particles  within  the  domain  so  that  their  centers  lie 
outside  the  volume  encompassed  by  other  particles.  Then,  forces  arising  due  to  particle  overlap  drives 
rearrangement  of  particle  positions  in  such  a  way  as  to  eliminate  overlap  while  continuously  settling 
under  the  influence  of  a  gravitational  force  [3].  Particles  are  discretized  using  constant  strain  triangu¬ 
lar  elements,  where  each  particle  contains  approximately  400-800  finite  elements  for  all  simulations 
performed  in  this  study,  depending  on  its  size.  The  combined  FEM/DEM  algorithm  has  been  verified, 
and  is  nominally  second-order  accurate  in  both  space  and  time. 


3  Predictions 

Preliminary  predictions  are  given  in  this  section  for  piston  supported  uniaxial  deformation  waves  in 
aluminized  HMX  having  (f>°s  =  0.835.  Macro-scale  wave  profiles  are  obtained  by  computing  spatial 
transverse  averages  over  meso-scale  fields  using  an  averaging  width  of  200  yum.  Though  computed 
average  profiles  do  vary  with  averaging  width,  numerical  experiments  indicate  that  meaningful  results 
are  obtained  provided  that  this  width  is  smaller  than  approximately  half  the  size  of  the  macro-scale 
wave  thickness. 

Figure  2  gives  predicted  variations  in  the  Hugoniot  curves  with  metal  mass  fraction  Anl. 
Figures  2(a)  and  (b)  express  these  curves  in  the  Up-PmiX  and  Up-D  planes,  respectively,  where  //„« 
is  the  average  equilibrium  mixture  pressure  behind  the  wave  and  D  is  the  corresponding  steady  wave 
speed.  To  this  end,  the  axial  dimension  of  the  computational  domain  was  chosen  sufficiently  long  to 
guarantee  that  steady  waves  evolved.  Here,  we  consider  piston  speeds  within  the  range  50  <  Up  < 
500  m/s  and  metal  mass  fractions  within  the  range  0  <  Am  <  1,  where  piston  speeds  in  excess  of 
approximately  Up  >  Up  =  250  m/s  result  in  nearly  complete  material  consolidation  for  all  values 
of  Am.  This  critical  value  of  piston  speed  separates  two  well-defined  impact  regions.  We  refer  to  the 
lower  pressure  region  corresponding  to  Up  <  Up  as  strength  dominated  because  material  strength  is 
important  in  preventing  complete  pore  collapse.  We  refer  to  the  higher  pressure  region  corresponding 
to  Up  >  Up  as  pressure  dominated  because  the  resulting  pressures  are  sufficient  to  completely  elim¬ 
inate  porosity.  Material  strength  plays  a  less  significant  role  than  pressure  in  the  pressure  dominated 
region. 

A  few  observations  are  noteworthy.  First,  Pmix  monotonically  increases  with  Am  due  to  the  high 
acoustic  impedance  of  aluminum  which  causes  the  reflection  of  acoustic  energy  at  explosive-metal 
interfaces  thereby  enhancing  the  deformation  and  pressure  of  explosive  particles.  This  observation 
is  commensurate  with  data  for  neat  granular  HMX  ( A„ ,  =  0)  and  neat  porous  aluminum  (A„,  =  1) 
[4].  Here,  neat  material  refers  to  one  that  consists  of  only  a  single  or  pure-phase  component.  Second, 
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Fig.  2.  Predicted  Hugoniot  curves  for  a  material  having  <f>°  =  0.835  and  metal  mass  fraction  Am:  (a)  the  P-U p 
plane;  (b)  the  D-Up  plane.  Here,  neat  refers  to  a  single  phase  granular  material. 


variations  in  D  with  Am  qualitatively  differ  between  the  two  impact  regions.  In  the  strength  dominated 
region,  D  slightly  decreases  with  increasing  Am  whereas  the  opposite  trend  is  predicted  for  the  pressure 
dominated  region.  Though  speculative  at  this  juncture,  this  prediction  is  plausible  because  the  lower 
aluminum  strength  enhances  dissipation  by  plastic  work  during  pore  collapse,  even  for  neat  granular 
aluminum  ( Am  =  1),  resulting  in  slightly  lower  wave  speeds  than  those  associated  with  neat  HMX 
in  the  strength  dominated  region.  However,  dissipation  by  plastic  flow  becomes  more  geometrically 
constrained  within  wave  profiles  in  the  pressure  dominated  region  due  to  the  elimination  of  porosity  re¬ 
sulting  in  higher  wave  speeds,  as  anticipated.  Though  not  shown  here  for  brevity,  the  weak  variation  in 
D  with  Up  in  the  strength  dominated  region  is  largely  due  to  inter-particle  friction  during  pore  collapse 
which  enhances  the  rigidity  of  the  material;  as  such,  we  also  refer  to  this  region  as  being  frictionally 
rigid.  Simulations  performed  for  frictionless  contact  indicate  a  larger  decrease  in  D  as  Up  — >  0  in  the 
strength  dominated  region. 

We  now  briefly  focus  on  predictions  for  wave  energetics.  Figures  3(a)  and  (b)  give  spatial  vari¬ 
ations  in  average  phase-specific  plastic  heating  rate  Wp  within  wave  profiles  corresponding  to  Up  = 
lOOm/s  and  Up  =  500  m/s,  respectively,  for  a  material  having  Am  =  0.2.  Here,  the  average  phase- 
specific  heating  rate  and  its  standard  deviation  are  computed  based  on  the  local  phase-specific  mass 
affected  by  plastic  work.  Also  shown  in  these  plots  for  comparison  is  the  equivalent  predicted  profile 
for  neat  HMX  having  the  identical  initial  configuration  as  the  aluminized  HMX  mixture.  Neat  HMX 
forms  a  proper  baseline  for  establishing  the  relative  impact  sensitivity  of  aluminized  HMX.  Average 
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Fig.  3.  Predicted  average  plastic  heating  rate  Wp  profiles  for  a  material  having  (f>°s  =  0.835  and  Am  =  0.2: 
(a)  Up  =  lOOm/s;  (b)  Up  =  500m/s.  Corresponding  standard  deviations  are  given  by  dashed  lines. 


heating  rates  within  wave  profiles  are  important  to  analyze  because  they  can  be  directly  compared 
to  predictions  given  by  macro-scale  multiphase  theories  and  because  they  give  a  good  indication  of 
deformation  wave  thickness. 

As  seen  in  Fig.  3(a)  for  Up  =  lOOm/s,  the  metal  plastic  heating  rate  exceeds  that  of  the  explosive 
within  the  wave  profile  for  aluminized  HMX.  Further,  the  explosive  plastic  heating  rate  is  marginally 
lower  than  that  of  the  equivalent  neat  HMX  configuration.  This  prediction  is  consistent  with  our  previ¬ 
ous  assertion  for  strength  dominated  waves  for  which  the  inclusion  of  aluminum  reduces  wave  speed 
by  enhancing  plastic  dissipation  within  aluminum  particles  during  pore  collapse,  thereby  reducing  ex¬ 
plosive  heating  relative  to  neat  HMX.  Standard  deviations  about  the  average  rates  are  large  due  to 
localized  plastic  deformation  near  inter-particle  contacts.  The  predicted  deformation  wave  thickness  is 
approximately  4  mm. 

Figure  3(b)  gives  predictions  for  Up  =  500  m/s.  This  case,  which  corresponds  to  a  pressure  domi¬ 
nated  wave,  results  in  a  marginally  higher  explosive  plastic  heating  rate  than  that  of  neat  HMX,  again 
possibly  due  to  enhanced  geometrical  confinement  associated  with  the  elimination  of  porosity.  The 
predicted  deformation  wave  thickness  is  approximately  1  mm.  The  implications  on  impact  sensitivity 
of  these  marginal  differences  in  explosive  heating  relative  to  that  of  neat  HMX  remain  an  ongoing 
topic  of  our  research. 

Figure  4  gives  corresponding  spatial  variations  in  average  phase-specific  friction  heating  rate  Wf 
within  wave  profiles.  A  friction  coefficient  of  /./  =  0.25  was  used  for  all  simulations  performed  in 
this  work  because  parametric  studies  indicate  that  this  value  results  in  maximum  frictional  dissipation. 
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Fig.  4.  Predicted  average  frictional  heating  rate  Wf  profiles  for  a  material  having  (p°s  =  0.835  and  Am  =  0.2: 
(a)  Up  =  lOOm/s;  (b)  Up  =  500m/s.  Corresponding  standard  deviations  are  given  by  dashed  lines. 


Here,  the  average  phase-specific  heating  rate  and  its  standard  deviation  are  computed  based  on  the 
local  phase-specific  mass  affected  by  friction  work,  which  is  approximated  by  finite-element  mass 
within  the  averaging  area  experiencing  sliding  contact.  Based  on  this  approximation,  friction  heating 
rates  exceed  plastic  heating  rates  in  all  cases  as  friction  work  is  a  surface  phenomenon.  Though  plastic 
work  affects  a  larger  mass  of  material,  friction  work  remains  a  significant  source  of  dissipation  for 
pressure  dominated  waves.  Importantly,  friction  work  is  believed  to  be  a  primary  mechanism  for  the 
formation  of  hot-spots  that  are  responsible  for  ignition  of  reactive  solids. 


4  Conclusion 

A  2-D,  Lagrangian,  combined  discrete  and  finite  element  technique  was  used  to  computationally 
characterize  piston  supported  uniaxial  waves  in  granular  metal-explosive  (aluminum-HMX)  mixtures. 
Predictions  indicate  that  the  impact  response  can  be  divided  into  material  strength  and  pressure  domi¬ 
nated  regions,  depending  on  piston  speed.  Strength  dominated,  low  pressure  waves  result  in  incomplete 
material  consolidation,  whereas  pressure  dominated  waves  are  sufficiently  strong  to  eliminate  porosity. 
Predicted  Hugoniots  indicate  that,  for  a  material  having  </>°  =  0.835,  the  wave  speed  slightly  decreases 
with  increasing  metal  mass  for  strength  dominated  waves  corresponding  to  Up  <  Up  =  250  m/s, 
whereas  the  opposite  trend  is  predicted  for  pressure  dominated  waves  corresponding  to  U p  >  Up. 
Inclusion  of  metal  is  found  to  marginally  decrease  the  average  explosive  heating  rate  within  strength 
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dominated  waves,  and  to  increase  explosive  heating  rate  within  pressure  dominated  waves.  Fric¬ 
tion  work  is  shown  to  be  a  potentially  significant  hot-spot  formation  mechanism,  even  for  pressure 
dominated  waves.  The  average  equilibrium  pressure  behind  waves  was  found  to  increase  with  metal 
mass  for  both  impact  regions  due  to  the  high  acoustic  impedance  of  the  metal.  It  is  emphasized  that 
careful  interpretation  of  plastic  and  frictional  heating  rates,  and  meso-scale  thermomechanical  fields, 
is  necessary  to  properly  characterize  the  thermal  response  of  reactive  solids. 
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